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Abstract 

Coupled 6-DOF/CFD trajectory predictions using an automated Carte- 
sian method are demonstrated by simulating a GBU-32/JDAM store sepa- 
rating from an F-18C aircraft. Numerical simulations are performed at two 
Mach numbers near the sonic speed, and compared with flight-test telemetry 
and photographic-derived data. Simulation results obtained with a sequential- 
static series of flow solutions are contrasted with results using a time-dependent 
flow solver. Both numerical methods show good agreement with the flight-test 
data through the first half of the simulations. The sequential-static and time- 
dependent methods diverge over the last half of the trajectory prediction, after 
the store produces peak angular rates. A cost comparison for the Cartesian 
method is included, in terms of absolute cost and relative to computing uncou- 
pled 6-DOF trajectories. A detailed description of the 6-DOF method, as well 
as a verification of its accuracy, is provided in an appendix. 
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1 Introduction 


Trajectory prediction is an important element in Computational Fluid Dynam- 
ics (CFD) simulations of bodies undergoing unconstrained, or partially constrained 
motion. Modeling this behavior involves integrating the Newton-Euler equations for 
six-degree-of-freedom (6-DOF) rigid-body motion, in response to aerodynamic and 
other externally applied loads. Numerous important applications for such models ex- 
ist, including store separation from an aircraft, booster separation from a space launch 
vehicle, canopy or shroud separation, and simulation of flight control systems. Many 
CFD technologies have been demonstrated for 6-DOF simulations, including struc- 
tured oversell, 2], unstructured tetrahedral[3, 4], and hybrid prismatic/Cartesian[5]. 
The current work demonstrates an integrated package for performing 6-DOF simula- 
tions couple with an inviscid, Cartesian embedded-boundary method. 

Such non-body-fitted, Cartesian methods are particularly interesting for 6-DOF 
applications since they can be made both extremely fast and robust, and the volume 
meshing can proceed automatically. Moreover, they are comparatively insensitive to 
the complexity of the input geometry since the surface description is decoupled from 
the volume mesh. In the current work, the “cut-cell” Cartesian meshing scheme of 
Aftosmis et al.[6] is utilized. The intersection of the solid geometry with the regular 
Cartesian hexahedra is computed, and polyhedral cells are formed which contain 
the embedded boundary. This volume meshing procedure is robust, computationally 
efficient, and does not require user intervention. 

In order to demonstrate the utility of the Cartesian 6-DOF package, a U.S. Navy 
GBU-32 Joint Direct Attack Munition (JDAM) store (cf. Fig. 1) separating from 
an F/A-18C is simulated using both sequential-static and time-dependent methods. 
This transonic JDAM separation -was put forward by the Navy as a “challenge” to the 
CFD community because it exhibited behavior that could not reliably be predicted 
with conventional store separation analysis tools (cf. Cenko [7, 8]) . The JDAM 
separation provides an attractive demonstration case because it contains a complex 
aircraft geometry, flight telemetry and photographic-derived quantitative data, and 
also because it has been simulated by numerous other CFD methods[9— 15], These 
previous CFD simulations can be broken into two broad classes; those which computed 
a set of static solutions which were used with a store trajectory simulation package, 
and those which computed the trajectory of the store within the CFD simulation 
process. Both of these approaches are supported with the current methods and a cost 
comparison will be presented. 

The discussion begins by reviewing the geometry used in the simulations, and 
briefly outlines the numerical scheme. Next it presents computed results for the 
JDAM separation flight conditions just below and just above sonic speed ( M ^ = 
0.962 and 1.055). These results are directly compared to both flight telemetry and 
photographic-derived data. The computational cost for the current method is pro- 
vided, along with a summary of the current results and topics for future work. A 
detailed description and verification of the stand-alone 6-DOF package used with the 
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Figure 1: U.S. Navy GBU-32 Joint Direct Attack Munition (JDAM) on the F/A-18C wing pylon. 
The dark green fins and center carriage provide the JDAM GPS guidance augmentation system 
which can be retrofit on a general purpose unit such as (in this case) an Mk-84. 


current scheme is included in an appendix. 


2 Numerical Scheme 

2.1 Geometry and Computational Mesh 

The surface geometry was provided as a set of structured surface patches. These 
were converted to water-tight surface triangulations of the various components. The 
addition of an internal duct connecting the engine diffuser face to the exit nozzle 
was required in order to form a water-tight fuselage. The component geometry for 
the complete F/A-18C is shown in Fig. 2, with water-tight components shown with 
different colors. All of the major components of the geometry are modeled, including 
the empennage, AIM-7 wingtip missile and rail, wing with leading-edge extensions 
(LEX), the LEX fence, the engine inlet including boundary layer vents, and the wing 
pylons holding a 330 gal. external fuel tank (EFT) inboard, and the GBU-32 JDAM 
outboard. Note that the flight configuration did not contain the AIM-7 wingtip 
missiles. Fig. 3 shows a closeup view of the JDAM in its initial position beneath the 
port, outboard wing pylon. The attachment hardware and ejector mechanism is not 

modeled. 

Using the automated Cartesian meshing scheme of Aftosmis et al.[6], the trian- 
gulated surface was used to generate an unstructured Cartesian volume mesh by 
subdividing the computational domain based upon the geometry. The sharp geo- 
metric features contain refined cells, while areas away from the geometry maintain a 
relatively coarse spacing. The intersection of the solid geometry with the the regular 
Cartesian hexahedra is computed, and polyhedral cells are formed which contain the 
embedded boundary. Regions interior to the solid geometry are removed. The solid- 
wall boundary conditions for the flow solver are then specified within these cut-cell 
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Figure 2: F/A-18C surface geometry. Water-tight components are shown with different colors. 



Figure 3: Closeup view of triangulated GBU-32/JDAM in its initial position beneath the wing 
pylon. 
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Figure 4: Cutting plane through the Cartesian volume mesh. 


polyhedra. 

In addition to mesh refinement near geometric features, pre-specified adaption 
regions are arranged around the major components of the F-18 aircraft to resolve the 
shock structures that occur at the current flow conditions. The adaption region which 
surrounds the JDAM translates with the center of mass (c.m.) location as the store 
drops. In the future, these pre-specified regions will be replaced with automated 
solution and geometry adaptation similar to the steady-state scheme outlined by 
Aftosmis and Berger[16]. A mesh refinement comparison was performed for the static, 
steady-state simulation with the JDAM in its initial position at the Mach 0.962 flight 
conditions (cf. Table 2). The resulting volume mesh is isotropic and contains 3.8M 
cells with a surface resolution of 1.0 in. A volume mesh cutting plane through the 
wing is shown in Fig. 4. Details of the mesh adaptation to the moving geometry will 

be presented in Sec. 3. 

2.2 Flow Solver 

The inviscid, parallel multigrid flow solver of Aftosmis et al. [17] provides static, 
steady-state flow simulations for Cartesian meshes. Recently, this flow solver has been 
extended to provide capability for time-dependent flows, including dynamic simula- 
tions with rigid bodies in relative motion[18, 19]. The current work implements an 
independent 6-DOF module which can be utilized as a stand-alone external applica- 
tion, or tightly coupled within the time-dependent flow solver. A flow diagram for 
the 6-DOF/CFD simulation process is shown in Fig. 5. The XML4CFD interface[20] 
is utilized to integrate the independent mesh generation, flow solver, post-processing, 
and 6-DOF steps into a unified computational framework. 
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New Geometry at n+1 


Figure 5: Process diagram for 6-DOF simulations. The red-colored processes are serial and the blue 
parallel. The XML4CFD interface[20] provides a single repository and API for the moving-body 
information required by the separate processes. 


The 6-DOF module decomposes the rigid-body motion into a translation of the 
center of mass and a rotation about an axis passing through the c.m. location. The 
position of the c.m. is updated using Newton’s laws of motion in the inertial frame, 
while the rotation of the body is determined by numerically integrating Euler’s equa- 
tions of motion in a body principal-axis system. The rotational position of the body 
is specified using Euler parameters, which are updated by numerical integration of 
the angular velocity. General external applied forces, in either the aerodynamic or 
body coordinate frames, can be specified. A detailed discussion of the 6-DOF model, 
along with validation test cases is presented in App. A. 

2.3 Ejector Force Model 

The JDAM is forced away from its wing pylon by means of identical piston ejectors 
located in the lateral plane of the store, -10.11 in. forward of the c.m., and 9.89 
in. aft. The ejectors extend during operation for 6 in., and the force of each ejector 
is a polynomial function of this stroke extension (cf. Cenko[7j). As the store moves 
away from the pylon it begins to pitch and yaw due to aerodynamic forces, and the 
stroke length of the individual ejectors responds asymmetrically. This response of 
the ejectors to the store motion is modeled, and the result is presented as a function 
of time for each piston. This modeling process for the F-18/JDAM is described by 
Fortin et al. [14] , and the results are presented in Table 1. A consistent theme with 
previous simulations is that this ejector model is limited (cf. [8, 9, 14, 15]). Since 
even slight errors in the initial trajectory of the store can become augmented as the 
separation simulation is marched forward in time, researchers have modified either the 
ejector model, or the computed JDAM trajectory, in order to provide a realistic store 
separation. Physically, the JDAM is constrained by the ejector mechanism, which is 
not accounted for in simplistic models. For example, the JDAM cannot be allowed 
to pitch nose-down without bound, as physically the aft ejector would restrict such 
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motions. 


Time(sec) 

Forward Ejector 
Force (lbf) 

Aft Ejector 
Force (lbf) 

0.00 

97 

97 

0.01 

206 

223 

0.02 

531 

283 

0.03 

1053 

549 

0.04 

4723 

988 

0.05 

4641 

4708 

0.06 

4542 

4633 

0.07 

4414 

4528 

0.08 

4255 

4386 

0.09 

0 

4243 

0.10 

0 

0 


Table 1: Modeled F/A-18C/GBU-32 ejector forces (cf. [7, 14]). 


While the focus of the current work is not to develop an ejector model for the F- 
18/JDAM configuration, simulating the store separation with an ejector model which 
has known inaccuracies serves little purpose. An attempt to modify the ejector model 
to account for the constraint imposed by the wing pylon and ejector mechanism is 
proposed. First, following the a posteriori observations of Cenko[8], the magnitude of 
the ejector forces was increased by 25%. Next, it is assumed that while the ejectors 
are accelerating (roughly 0.0 < t < 0.05 sec), the rotation of the JDAM is restricted 
by friction between the ejector pistons and the JDAM surface. From examining the 
flight data it is clear that the rotation is not completely restrained, so a friction 
resistance equivalent to 50% of the aerodynamic moments is imposed initially, which 
is allowed to linearly decrease to no resistance at t. = 0.05 sec. This modified ejector 
model is used with all the simulation results presented here. 


3 Computed Results 

The numerical scheme outlined in the previous section was used to compute the 
separation of a GBU-32/JDAM from an F/A-18C at the two flight conditions listed 
in Table 2. The inertial properties for the JDAM were provided by the Navy, and 
are summarized in Table 3. The pylon ejector modeling was discussed in Sec. 2.3. 
This configuration was tested in the wind-tunnel using a Captive Trajectory System 
(CTS) and in-flight by the U.S. Navy (cf. Cenko[7]). Near sonic speeds, the variation 
of pitching and yawing moments experienced by the JDAM with Mach number be- 
comes highly non-linear. This strong non-linearity makes trajectory prediction using- 
linearized methods (cf. Keen[21]) challenging. High-fidelity CFD methods can po- 
tentially provide a cost-effective, accurate tool for predicting store trajectories at all 
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flight conditions. 



Case 1 

Case 2 

Mach number (Moo) 

0.962 

1.055 

Altitude (h) 

6332 ft. 

10,832 ft. 

AOA (a) 

0.46° 

-0.65° 

Dive Angle (7) 

43.0° 

44.0° 


Table 2: Computed flight conditions. 


S re f 

1.767 sq. ft. 

■^ref 

1.5 ft. 

c.m. 

62.66 in. from nose 

mass 

2059.44 lbm j 

1 XX 

20.02 slug - sq. ft. 

lyy 

406.56 slug - sq. ft. 

h. 

406.59 slug - sq. ft. 

hz 

-0.680 slug - sq. ft. 

Ixy 

0.860 slug - sq. ft. 

_k_j 

0.00 slug - sq. ft. 


Table 3: GBU-32 JDAM inertial properties and reference quantities. 

Static, steady-state simulations were computed with the JDAM in its initial po- 
sition below r the wing pylon for both flight conditions. Surface pressure contours on 
the body surface are showrn in Fig. 6 for the M <*, = 1.055 simulation. The shock 
reflections on the wing pylons due to the stores are visible, as are the shocks that 
appear on the canopy, wing, and empennage. The cutting plane shows the resolution 
of the shocks to the farfield. 

The computed forces and moments on the JDAM from the initial static simu- 
lations are compared with wind-tunnel and flight data in Table 4. No uncertainty 
predictions or error estimates are available for the wind tunnel or flight data. The 
computed results are in good agreement with the flight and tunnel data, with the 
largest discrepancy occurring in yawing moment at M 0 0 — 0.962, which is less than 
10% variation. In general, the computed results compare more favorably to the flight 
data at M x = 1.055 than 0.962, as would intuitively be expected. 

3.1 Sequential-static Simulations 

The current wwk simulates the separation of the JDAM using both time-dependent 
and steady-state methods. The inertia of the GBU-32 is very large, and the expecta- 
tion is that unsteady effects are minimal, at least while the store is still close to the 
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Figure 6: Surface pressure 


contours on the F-18C surface (Moo — 1 .055, a — 0.65 ). 


pylon. This thesis is examined by comparison of time-dependent separation results 
with ’’sequential-static” simulations. The sequential-static results are presented first. 
In this method, the store is repositioned at the new time level based upon the com- 
puted loads at the previous time level (cf. Fig. 5, Sec. 2.2), however the flow solver 
ignores the motion of the body and treats it as a static, steady-state problem at the 
new body position. This approach can be attractive when accurate, time-dependent, 




1 

c A 

c Y 

C N 

Ci 

c m 

c n 

Wind Tunnel 

- 

0.31 

0.11 

- 

-2.32 

-2.76 

Flight 

- 

1 0.31 

0.15 

- 

-2.5 

-2.8 

Computed 

0.67 

0.33 

0.09 

0.16 

-2.36 

-2.49 


a) Moo = 0.962, a = 0.46°, 7 = 43° 



C.4 

Cy- 

Ov 

Ci 

c 

c n 

Wind Tunnel 


0.24 

-0.02 

~ 

-2.07 

-2.56 

Flight 

- 

0.25 

-0.05 

~ 

-2.0 

-2.2 

Computed 

0.65 

0.28 

-0.03 

0.15 

-2.02 

-2.11 


b) Moo = 1 055, a = -0.65°, 7 = 44° 


Table 4: Computed forces and moments on the JDAM for the initial store position. Wind tunnel 
and flight data taken from Cenko[7], 


moving-body flow solvers are not available. In the current work, the computed so- 
lution at the previous time level is transfered to the new mesh, after the body has 
been repositioned, to use as an initial guess. This transfer process, which is described 
in [19] for the time-dependent scheme, minimizes the computational cost since the 
solution at the previous time level provides a good initial guess for the solution at 
the next time level. 

A constant timestep of At = 0.0075 sec. is used for these simulations. Due to 
time constraints, it was not possible to perform a time resolution study for these 
cases. Information travels roughly one JDAM body length in 12 timesteps using this 
resolution, which is felt to be reasonable. All simulations were run through t = 0.45 

sec. 

Computed results for the relative displacement of the JDAM c.m. location are 
compared to flight data for both computed cases in Fig. 7. Similar plots for the 
angular position and angular velocity of the JDAM are shown in Figs. 8 and 9 re- 
spectively. Below t = 0.20 sec. the predicted displacement and angular position are 
in good agreement with the flight data, however the angular rate prediction has be- 
gun to degrade. At later times, the cumulative errors in angular position lead to a 
poorer agreement with the flight data, while the predicted displacement of the c.m. 
correlates well through the simulation. The accuracy of the current predictions is 
commensurate with previous computed results for this same configuration[9-15]. The 
degradation of the predicted angular orientation will be discussed in the next section 
with the time-dependent simulations results. 

The miss distance, or the distance between the closest points on the JDAM and 
any other component of the aircraft, is presented in Fig. 10. While the predicted 
displacement and angular position are in good agreement with the flight data over 
the time interval presented, the miss distance underpredicts the separation between 
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(a) Mac = 0.962, a = 0.46°, 7 = 43° (b) M x = 1.055, or = -0.65°, 7 - 44° 

Figure 7: Relative displacement for sequential-static simulations. Flight data from Cenko[7], 



(a) Moo = 0.962, a = 0.46°, 7 = 43° (b) M x = 1.055, a = -0.65°, 7 - 44° 

Figure 8: Angular positions for sequential-static simulations. Flight data from Cenko[7], 

the store and the wing pylon. The explanation for this is that as the ejectors push 
away the store, there is a reaction force applied to the pylon. This reaction leads 
to a rolling moment on the aircraft which rolls the pylon away from the JDAM, i.e. 
increases the miss distance between the two. This reaction of the aii craft is not 
modeled in the current work (or in previous work in the literature), and hence the 
separation is underpredicted. The closest miss, which occurs near t = 0.10 sec., is 
caused by the tail fins sweeping under the pylon as the JDAM yaws nose outboard. 
At t — 0.20 sec. the closest component changes from the pylon to the EFT, as the 
body continues to yaw and fall. 

Figure 11 shows a series of snapshots of the surface pressure as the JDAM falls 
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(a) Moo = 0.962, a = 0.46° , 7 = 43° 


(b) Moo = 1 055, q = -0.65°, 7 = 44° 


Figure 9: Angular rates for sequential-static simulations. Flight data from Cenko[7]. 




(a) Moo = 0.962, a = 0.46°, 7 = 43° (b) = 1.055, a = -0.65°, 7 = 44° 


Figure 10: Miss distances for sequential-static simulations. Flight data from Cenko[7], 

through t = 0.40 sec. in the M* = 1.055 simulation. The nose of the store is forced 
downward and outboard by the shock from the leading-edge of the wing. This causes 
the JDAM to pitch and yaw immediately upon release from the holding pylon. The 
change in shock structure on the pylon as the JDAM releases can be seen, as well 
as changes on the aft portion of the aircraft fuselage. As the JDAM falls, the tail 
fins provide restoring moments which cause the store pitch back nose up and inboard 
(compare with Figs. 8 and 9). A complementary series of snapshots which show the 
adaptation of the mesh to the moving geometry are shown in Fig. 12. The mesh 
automatically adapts to the new geometry position, and also coarsens in regions the 
body has moved through. 
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(a) t = 0.0 sec. 


(b) t — 0.1 sec. 




(c) t = 0.2 sec. 


(d) t = 0.3 sec. 


Figure 11: Surface pressure contours during JDAM separation. 


3.2 Time-dependent Simulations 

The previous section presented results of coupled 6-DOF/CFD trajectory pre- 
dictions using sequential-static flow simulations. This is contrasted here with fully- 
coupled, time-dependent trajectory simulations performed using the Cartesian moving- 
body solver described in [19]. Analyzing Fig. 9, the angular rate prediction for the 
sequential-static simulations begins to degrade after the rotation of the body experi- 
ences both the highest velocities and an inflection point in the acceleration, i.e. near 
t = 0.125 for pitch rate, and t = 0.20 for yaw rate. This combination of high velocity 
and change in sign of acceleration indicate regions in the store trajectory where dy- 
namic, or unsteady effects, may be significant. This is examined in Figs. 13-16, which 
present relative displacement, angular orientation, angular rate, and miss distance for 
the time-dependent simulations, compared with the sequential-static simulations and 
flight-test data. The data shows that the two CFD trajectory simulations are in good 
agreement prior to t = 0.125, when the pitch rate reaches a maximum. After this 
point, the predicted pitch behavior is improved, however the yaw prediction degrades. 
The pitch and yaw trajectories are similar in the sequential-static and time-dependent 
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(a) t — 0.0 sec. 


(b) t — 0.1 sec. 



(c) t = 0.2 sec. 


(d) t = 0.3 sec. 


Figure 12: Cutting pianes through the volume mesh during JDAM separation. 


simulations, except for the response near the maximum rates, i.e. the dynamic ef- 
fects are largely localized to this region of the trajectory. The relative displacement 
prediction is nearly unchanged in the time-dependent simulations at = 1.055, 
however M x = 0.962 shows a relatively significant change in vertical drop, which is 
not currently well understood. The underprediction of the separation distance after 
t — 0.20 is caused by the over-predicted yaw angle in both the sequential-static and 
time-dependent simulations, which causes the tail fins to remain close to the EFT. 

Consistently, in both the sequential-static and time-dependent simulations, the 
predicted roll behavior of the JDAM does not correlate well with the flight data. 
This is not unique to the current work, and has been noted in previous trajectory 
predictions for this configuration [7 -15]. Cenko[7] notes “[roll attitude] is the hardest 
to predict, fortunately has a minimal impact on the trajectory”. While it’s true 
that small changes in roll orientation are likely insignificant, the current predictions 
consistently vary from the flight data by roughly 5° of roll, and even while the store is 
still being pushed by the ejectors the roll is predicted in the opposite direction. Since 
the roll orientation can effect the restoring moment provided by the tail fins, it’s 
unclear whether these small differences can accumulate to produce the larger errors 
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(a) Moo = 0.962, a = 0.46°, 7 = 43° 


(b) Moo = 1 055, a = -0.65°, 7 = 44° 


Figure 13: Relative displacement for time-dependent simulations. Flight data from Cenko[7], 



Figure 14: Angular positions for time-dependent simulations. Flight data from Cenko[7], 


seen in pitch and yaw prediction in the current work. 

3.3 Computational Cost 

The computational cost for the current Cartesian/6-DOF scheme is presented in 
two forms; absolute and relative to computing a fixed database of static results. 
Note that the current work was performed with tools designed for computing a single 
fixed static simulation, and little effort has gone into tailoring them for sequential 
moving-body calculations. All simulated results presented here were computed using 
NASA Ames’ 1024 CPU, single-image SGI Origin 3000 (03K) which has 600MHz 
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(a) Moo = 0.962, a = 0.46°, 7 = 43° (b) M x = 1.055, a = -0.65°, 7 = 44° 


Figure 15: Angular rates for time-dependent simulations. Flight data from Cenko[ 7 ]. 




(a) Moo = 0.962, a = 0.46°, 7 = 43° (b) M x = 1.055, a = -0.65°, 7 = 44° 


Figure 16: Miss distances for time-dependent simulations. Flight data from Cenko[ 7 ], 


MIPS4 processors. The current flow solver has been demonstrated to scale linearly 
to 512 CPUs on this architecture for problems of the size considered here. The 
current simulations all required roughly 260 single-CPU-hours of computational time 
to complete, with less than 5% of the computational time utilized by the volume mesh 
generation process. The sequential-static and time-dependent simulations require the 
same computational time with the current scheme. The wallclock time to complete a 
simulation using 32 CPUs is approximately 15 hours. This time reflects the adverse 
effects of the serial mesh generation on the parallel efficiency of the entire process. 
Parallelizing the entire process, including volume mesh generation, will be a major 
focus of future work. 
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The current work couples together the CFD flow solver and the 6-DOF trajectory 
prediction. Another method of integrating high-fidelity CFD with 6-DOF predictions 
is to build a computational database of results, and then “fly’ 6-DOF trajectories 
through this computed database. The advantage of this approach is that once the 
initial database is created, many 6-DOF trajectories can be computed essentially at 
no cost. The disadvantage of the database approach for 6-DOF simulations is the 
large number of computational cases required to build even a minimal database. For 
a single fixed wind vector (Mo c , a, ( 3 ) there are 6 free parameters (3 displacement 
and 3 angular position) for a static CFD database.* If each of these is allowed to 
vary over 10 distinct states (which is relatively coarse), then 10 6 computed cases are 
required to fill the database. This is impractical even for wind tunnel programs. It’s 
possible to reduce the required independent variables by assuming that the horizon- 
tal and lateral relative displacements are much less than the vertical, and that the 
roll orientation of the body can be ignored. This reduction leaves on the order of 
1000 data points required for steady-state simulation. In the current work, an ini- 
tial steady-state calculation is required at the initial position of the store, and each 
timestep costs roughly 1/5 of a full static simulation. As 60 steps were required for 
a full simulation using the current timestep, the cost for the current coupled 6-DOF 
trajectory simulations is roughly 10 steady, static simulations. This implies that on 
the order of 100 such coupled simulations can be performed for the cost of building 
a coarse, approximate database. Further, each coupled simulation is independent, so 
that the simulations can be carried out in parallel. The higher accuracy and relatively 
low cost makes these coupled CFD/6-DOF simulations an attractive analy sis tool. 


4 Summary 

The utility of a coupled Cartesian/6-DOF trajectory prediction scheme has been 
demonstrated by simulation of a GBU-32 JDAM separating from an F/A-18C. The 
Cartesian scheme provides an automated, robust meshing scheme which can easily 
be integrated into a design analysis. The accuracy and computational cost of the 
current simulated results are commensurate with previous results for the F-18/JDAM 
separation computed using body-fitted approaches. 

Future work will progress on two major fronts; understanding the discrepancies 
in predicted angular orientation that occur at later time levels, and optimizing the 
flow simulation process for these moving-body simulations. There are many possi- 
ble explanations, both computational and experimental, for the degradation in the 
predicted trajectory at later time levels. It’s important to understand whether this 
behavior is related to the current approach so that it can be corrected, if necessary. 
The process optimization itself will mainly focus on parallelizing the volume mesh 
generation, and incorporating an solution-adaptive capability. 

‘Static here refers to the absence of any dynamic stability derivative information. 
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Appendix 


A 6-DOF Model 

This appendix describes an implementation of the unconstrained motion of a 
rigid body, commonly referred to as six-degree-of-freedom (6-DOF) motion. The 6- 
DOF model is implemented as a stand-alone package with a well-defined Application 
Programming Interface (API). In this manner it can easily be integrated within a CFD 
flow solver, or similar application, or used as a stand-alone package, for example when 
performing trajectory simulations within a pre-existing database of force and moment 
data,. 

The 6-DOF motion is computed by solving the Newton-Euler equations for rigid- 
body motion. The motion is broken into a translation of the center of mass (c.ra.) 
of the body (Newton’s equations), and a rotation about a centroidal axis system 
attached to the body (Euler’s equations) (cf. Fig.A.l). Here superscripts are used to 
designate the coordinate system, with i referring to the inertial frame, and b the body 
frame. The inertial frame is considered to be the natural coordinate system of the 
geometry. Note that this inertial frame is not in general identical to the aerodynamic 
frame in which forces and moments are calculated, so a transformation from the 
aerodynamic frame to the inertial frame is required. 



Figure A.l: Inertial and body-fixed coordinate systems. Superscripts are used to designate the 
coordinate system, with i referring to the inertial frame, and 6 the body frame. The body frame 
is rotated by an angle <b about the axis a relative to the inertial frame. The body frame is the 
unique frame defined by the principal axes of the moments of inertia. This is contrasted with the 
non-unique body frame b' which is defined by convenience. The angular velocity is u> is the principal 
axes frame, and (p, q , r) in the general body frame. 

The mass center translation is governed by Newton’s laws of motion, which are 
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written in the inertial frame as 

F = Fj + F‘ + F; = mr’. n . ( A1 ) 

where the applied force acting through the center of mass has been broken into three 
components; the aerodynamic forces F^, the external applied forces (such as thrust) 
F* , and the forces due to gravity F^. Equation A.l is written in non-dimensional 
variables using the reference density (p a Q ), reference velocity (freestream some speed 
a oo ) j and a reference length ( L ). The non-dimensional mass is thus the dimensional 
mass scaled by the mass contained in a reference unit volume 

til 

Poo L* 

and similary the forces and gravity are non-dimensionalized by 

F 

F P^a^L 2 


Newton’s laws can be integrated directly to give the position of the mass center as 
a function of time. Holding F constant over the discrete physical timestep (t , t ) 
gives . 

4m (t n+1 ) = l — Ai2 + U C .m.(t n )At + 4 m (t") (A.2) 

cm - 2 m 

where u). m is the velocity of the center of mass. 

The rotational motion is governed by Euler’s equations of motion. The body axes 
are specified to coincide with the principal axes of inertia, with origin at the center 
of mass (cf. Fig. A.l). Euler's equations are then 


Mi = I b A - (I b 2 - / 3 V 2^3 

M b = I b 2 dj b - (I b 3 - I b M ( A - 3 ) 

M b = I b u b -(I b -I b M4 


where M 6 are the applied moments in the body frame, and are broken into aero- 
dynamic and external components as in Eqn.A.l. u> b is the angular velocity in the 
body frame, and I 6 are the principal moments of inertia. Using the same reference 
quantities as above, the non-dimensional applied moments and moments of inertia 
are given by 


M 6 

I 6 


M fe 

PocO^L 3 

I b 

~pj> 
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Equation A. 3 is integrated numerically using a 4 t/l -order Runge-Kutta scheme. 

In order to transform the angular velocity into a change in orientation, it’s desir- 
able to use quaternions, commonly referred to as “Euler parameters” , to specify the 
angular orientation of the body frame with respect to the inertial frame (cf. Fig. A.l). 
A quadratic transformation matrix A, which is composed of the direction cosines, can 
be expressed as the result of two successive linear transformations 

A = GL r 

where G and L 7 are both composed of the 4 Euler parameters 

p = [e 0 ei e 2 e 3 ] T 

The transformation matrix in terms of the Euler parameters is given by 

e o + e ? - k e i e 2 - e 0 e 3 eie 3 + e 0 e 2 ‘ 

A = 2 eie 2 + eoe 3 e o T | e 2 e 3 — eoei (A. 4) 

Cie 3 — eo e 2 e 2 e 3 + eo e i T c 3 — \ 

The Euler parameters specify an axis of rotation (a), and an angular displacement 
about that axis (d>) 

d> 

e 0 = cos - 

• <P 

ei - a x sin — 

l < A5 ) 

e 2 = CLy sin - 
y 2 

. 0 

e 3 = a z sin - 

According to Euler’s theory of motion, the Euler parameters are the same in both 
the body and fixed reference frames, so no superscript appears on p, however note 
that in this case the discussion assumes the reference frame is attached to the center 
of mass. 

Using this, the change in orientation due to rotation can be found through 

P = ^L r u, 6 (A. 6) 

which can also be integrated numerically using a 4 t/l -order Runge-Kutta scheme. Since 
the Euler parameters are unit-normalized quaternions, it’s necessary to impose that 
|p| = 1 after solving Eqn. A. 6. L T is given by 

~ e l ~^2 ~ e 3 

jT _ e o — e 3 e 2 

e 3 e 0 - &\ 

— e 2 Ci Cq 
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In order to update the position of an unconstrained rigid body, the following 
procedure is thus followed 

1. Translate 

Solve Newton’s law’s of motion using Eqn. A. 2 for the translation of the center 
of mass. 

2. Rotate 


(a) Angular Velocity 

Numerically solve Euler’s equations of motion using Eqn. A. 3 for the an- 
gular velocity in the body frame. 

(b) Euler Parameters 

Update the orientation of the body by numerically integrating Eqn. A. 6 
for the Euler parameters. 


3. Reposition 

Position the body according r‘ = r’ m 4- Ar b 


While quaternions are convenient for calculating the angular position of a rigid 
body, they are not always intuitive. It’s often desirable to transform the quaternions 
to a set of three angles which are non-unique, but in practice often unambiguous. 
The Euler parameters can be converted to angular displacement ((f) x , (f> y . (f> z ) of the 
body relative to the inertial frame using 


2.0 (eoei + 6263) 

tan(0 x J — 2 2 


en — e? — ei, + 


-0 c i 

sin (4> y ) = — 2 . 0 (eie 3 — e<je 2 ) 

; 2.0 (eie2 + eoes) 

tan(©.J = e2 , e 2 _ e 2 _ e 2 


(A.7) 


The potential singularity in x and z orientation is obvious. 


A.l 6-DOF Model Verification 

The 6-DOF implementation was verified using a variety of analytic test cases. 
Integration of Newton’s laws (Eqn. A.l) is verified using response of a point mass to 
a constant external force, and the terminal velocity of a falling sphere examines this 
integration for a non-constant external force. Integration of Euler’s Eqn. A. 3 was 
examined using the response of a cylinder undergoing a coupled spin. A tumbling 
rectangular volume demonstrates that the numerical implementation has the same 
stability properties as the physical system. 

Translation is integrated analytically according to Eqn. A.l, holding the applied 
force constant over the timestep. Since the integration is analytic, it is exact in the 
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presence of a constant external force. To verify this consider a point mass with an 
initial upward velocity in a gravitational field. The force due to gravity is normalized 
such that mg = (0,0, —1), and the initial velocity vector is u = (0,0, 1). Figure A. 2 
shows the exact solution for this case compared with the computed solutions taken 
with At = 0.025 and At. = 0.1. Since the integration is exact for this constant external 
force, the numerical integration reproduces the exact solution at both timesteps. 



Figure A. 2: Distance as a function of time for a unit point mass in a gravitational field. The force 
due to gravity is normalized such that mg = (0, 0, —1), and the initial velocity vector is u = (0, 0, 1). 

When the external force is non-constant, holding F fixed over the timestep results 
in formal first-order accuracy. To demonstrate this, consider a sphere with drag 
coefficient Cp = 0.5 falling through air in a gravitational field. The external force is 
F = mg — Cp\pu 2 z S. When the gravity and drag forces balance each other the sphere 

reaches its terminal velocity «<*, = • Taking p = 1 , S = 1 , and mg = (0, 0, —1) 

to construct a unit model problem, one can solve for the velocity as a function of 
position for an object initially at rest and falling in the —2 direction as 

u z {z) = v/4(l -e 0 - 5 *) 

Figure A. 3 plots velocity as a function of distance for the theoretical result and 
numerical experiments run with At = 0.1, 0.2 and 0.4. As expected halving the 
timestep halves the maximum error in the simulations providing the expected order 
of accuracy. All simulations converge to the correct terminal velocity (u ^ = —2.0) 
since the external force becomes a constant at this limit. 

Smart [1] presents an exact solution to Euler’s equations of motion that corre- 
sponds to a tumbling and spinning cylinder. For this example there is no translation, 
and in the current notation the inertial properties are given by 

h = h 

I\ — h = h ~ h = a h 
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Figure A. 3: Velocity vs distance evolution of a sphere falling in a gravitational field subject to air 
resistance. 

The analytic solution of the time evolution of the angular velocities is 

uj 1 = ucos(Af) 
u >2 = 6sin(Af) 

u>3 = C 

Setting Ii = h = 1 . 0 , ct = 0.5, a = 1.0 and c = 0.5 gives / 3 = 0.5 and A = 0.25. 

At t = 0, the 6-DOF model is initialized with w = (1.0, 0.0, 0.5). Figure A.4 
shows the system’s response to this initial condition compared against the analytic 
solution with b = — 1 . Evolution of the numerically-integrated angulai velocities are 
plotted for At = 1.0, 2.0 and 4.0. Since the integration is formally fourth order, 
the results convei'ge very quickly and only the symbols for At = 4.0 clearly diffei 
from the theoretical curve. With a At of 1.0 there are approximately 12 samples pei 
wavelength. Table A.l provides a quantitative comparison of the convergence, listing 
the error in uji at t = 100 as At increases. Since the error increases with time, this is 
the maximum error or the interval t = [0, 100]. The data in Table A.l shows fourth 
order asymptotic convergence, as expected. 
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Figure A. 4: Time evolution of numerical integration for angular velocities compared with an exact 
solution of Euler’s laws of motion from Smart [1]. 


At 

%Error in u\ 
at t — 100 

Improvement 

Order of accuracy 

0.25 

5.89E-5 

20.96 

4.39 

0.50 

1.2047E-3 

22.80 

4.51 

1.0 

0.02747 

25.08 

4.65 

2.0 

0.6890 

24.09 

4.59 

4.0 

16.60 

- 



Table A.l: Accuracy of numerical integration of Euler’s Eqs. of motion for coupled rotation of a 
cylinder. 


While the system of Euler’s Eqns. decouples when the rotation axes are aligned 
with any one of the principal axes of a body, stability analysis shows that this rotation 
is only stable around the minor or major axis - rotation around the semi-major axis 
is unstable. The coupling of the system means that any small perturbation about 
the semi-major axis will excite rotation about the others (cf. Thompson[2]). With 
I\ = 1, 1-2 = 10, and I 3 = 100, the 6-DOF model was initialized with u> = (0,0, 1), 
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prescribing rotation around the major axis. The system was perturbed by imposing a 
moment with magnitude 0.01 about the minor axis over the first time step (At — 0.1). 
Figure A. 5 shows the system response in terms of the angular rates around the minor 
and semi-major axis. As expected, the initial perturbation excites oscillations around 
both of these axes, but these oscillations disappear rapidly as the system stabilizes. 
Since the system is lossless (i.e. contains no physical dissipation), the rotational energy 
of the system must be conserved. Figure A. 6 shows the Euler angles of the object, 
revealing that the initial oscillations are transformed into a steady, but extremely 
small, oscillation about both the minor and semi-major axes. Figure A. 7 shows that 
this oscillation persists undamped, as expected from a lossless system. 



Figure A. 5: Time evolution of angular velocity around minor (red) and semi-major (blue) axes for 
a system spinning around major axis. System is perturbed at t = 0 with an impulsive couple around 
the minor axis. 

Contrast the results of Figs. A.5-A.7 with those shown in Fig. A. 8. In the example 
shown in Fig. A. 8, the initial angular velocity is prescribed as w = (0, 1,0) and the 
moments of inertia are unchanged from the previous example. Spin is therefore around 
the semi-major axis. When the same initial perturbation is applied, the perturbation 
is amplified, resulting in spin around all three axes. It's clear from the plot that the 
magnitude of the resulting angular rates is proportional to the moments of inertia 
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Time 


Figure A. 6: Time evolution of Euler angles showing small oscillation excited by perturbation of 
system spinning around major axis. System is perturbed at t — 0 with an impulsive couple around 
the minor axis. 

around these axes. Again, since the system is lossless, this tumbling behavior persists 
undamped. 
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Figure A. 8: Angular rate response to an initial perturbation for an object initially spinning around 
its semi-major axis at rate UJ 2 — 1. Coupling quickly leads to a tumbling motion. 
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